/****************************************************************************************************
* Title: Workhorses of Opportunity, Howard and Weinstein
*
* Inputs:
*   - justnormasylum.dta
*   - county_outcomes_dta.dta
*
* Outputs:
*   - xxx_with_average_by_race.pdf (one for each outcome variable)
*
* Description:
*   This script estimates the effect of being from a "normal school" county on various outcomes
*   across racial groups and parental income percentiles. It generates comparison graphs of these 
*   effects with confidence intervals and baseline averages.
****************************************************************************************************/

* Load and prepare normal school indicator
clear
use "justnormasylum"
rename hasnormalschool hasnormal
keep hasnormal cty_fips 
tempfile normal
save `normal'

* Set up environment
clear
clear matrix
clear mata
set maxvar 15000

* Load Chetty county-level outcome data and merge with normal school indicator
clear
use county_outcomes_dta
gen cty_fips = state*1000 + county
rename state statefips
rename county countychetty
merge 1:m cty_fips using `normal', keep(3)

* Set up outcome and result placeholders
gen outcome = .
gen outcome0 = .
tempfile results

* Loop over variables, percentiles, and racial groups to estimate treatment effects
foreach xxxxx in has_dad has_mom hs coll comcoll somecoll married jail staycz stayhome working kfr kir {
    foreach percentile in 1 25 50 75 100 {
        foreach race in white black hisp asian {
            
            replace outcome = `xxxxx'_`race'_pooled_p`percentile'
            replace outcome = outcome * 100 if inlist("`xxxxx'", "kfr", "kir")
            
            reghdfe outcome 1.hasnormal, absorb(statefips) cluster(statefips) nocon
            
            local numobs = e(N)
            summ outcome if !hasnormal
            local basemean = r(mean)
            
            preserve
                parmest, fast
                gen percentile = `percentile'
                gen race = "`race'"
                gen outcome = "`xxxxx'"
                gen baselinemean = `basemean'
                gen numobs = `numobs'
                cap append using `results'
                save `results', replace
            restore
        }
    }
}

* Load and process results
clear
use `results'
gen min90 = -invttail(dof, .05)*stderr + estimate
gen max90 = invttail(dof, .05)*stderr + estimate

* Create and export race-stratified graphs for each outcome
foreach xxxxx in hs coll somecoll married jail staycz stayhome working kfr kir {
    preserve
        keep if outcome == "`xxxxx'"
        
        * Offset for jittering by race
        replace percentile = percentile - 1 if race == "white"
        replace percentile = percentile + 1 if race == "hisp"
        replace percentile = percentile + 2 if race == "asian"
        
        twoway scatter estimate percentile if race == "white", mcolor(dkgreen) yline(0, lpattern(dash) lcolor(gs8)) ytitle("Effect of Normal County (Spikes)") xtitle("Parents' Income Percentile") xlabel(0 25 50 75 100) || /// 
            rspike min95 max95 percentile if race == "black", color(maroon) || ///
            rcap min90 max90 percentile if race == "black", color(maroon) || ///
            scatter estimate percentile if race == "black", mcolor(maroon) || ///
            rspike min95 max95 percentile if race == "hisp", color(navy) || ///
            rcap min90 max90 percentile if race == "hisp", color(navy) || ///
            scatter estimate percentile if race == "hisp", mcolor(navy) || ///
            rspike min95 max95 percentile if race == "asian", color(orange_red) || ///
            rcap min90 max90 percentile if race == "asian", color(orange_red) || ///
            scatter estimate percentile if race == "asian", mcolor(orange_red) || ///
            line baselinemean percentile if race == "white", lcolor(dkgreen*.5) lpattern(dash) yaxis(2) ytitle("Average Value in Asylum County (Dashed)", axis(2)) name(`xxxxx'_average, replace) || ///
            line baselinemean percentile if race == "black", yaxis(2) lcolor(maroon*.5) lpattern(dash) || ///
            line baselinemean percentile if race == "hisp", yaxis(2) lcolor(navy*.5) lpattern(dash) || ///
            line baselinemean percentile if race == "asian", yaxis(2) lcolor(orange_red*.5) lpattern(dash) || ///
            rspike min95 max95 percentile if race == "white", color(dkgreen) || ///
            rcap min90 max90 percentile if race == "white", color(dkgreen) ///
            legend(order(1 4 7 10) label(1 "White") label(4 "Black") label(7 "Hispanic") label(10 "Asian"))
        
        graph export "`xxxxx'_with_average_by_race.pdf", as(pdf) replace
    restore
}
